/**
 *  @file
 *  @brief   Hexa8SRI单元
 *  @author  xiashijun
 *  @date    2025.03.11
 *  @version v1.0
 *  @par Copyright(c):  ShenGong  Simulation corporation
 */
#include "Hexa8SRIElementCalculator.h"

#include "Utility/Logging/SGLogger.h"

#include "DataStructure/Input/Element/CHEXAData.h"
#include "DataStructure/Input/Solver/DOFDirections.h"
#include "DataStructure/Input/Solver/TenserTranform.h"

#include "DBManager/IService/Input/Element/IElementService.h"
#include "DBManager/IService/Input/Material/IMaterialSummaryService.h"
#include "DBManager/IService/Input/Node/INodeService.h"
#include "DBManager/IService/Input/Property/IPSOLIDService.h"
#include "DBManager/IService/Result/ElementalStrain/IHEXAStrainService.h"
#include "DBManager/IService/Result/ElementalStress/IHEXAStressService.h"

#include "SGFem/Pre/Assembly/Assembly.h"
// #include "SGFem/Pre/Assembly/CoordinateProcessor.h"
#include "SGFem/Pre/Assembly/ModelInfoQuery.h"

#include "../SolidElement/SolidMaterialParser.h"

using SG::Algebra::Matrixd;
using SG::Algebra::PointArray;
using SG::Algebra::Vector3D;
using SG::DataStructure::Common::Id;
using SG::DataStructure::Common::IdArray;
using SG::DataStructure::Common::Real;
using SG::DataStructure::FEM::CHEXA8Data;
using SG::DataStructure::FEM::DOFDirections;
using SG::DataStructure::FEM::DofSet;
using SG::DataStructure::FEM::ElementTypes;
using SG::DataStructure::FEM::PSOLIDData;
using SG::DataStructure::FEM::TenserTranform;
using SG::DataStructure::Post::HEXAStrainData;
using SG::DataStructure::Post::HEXAStressData;
using namespace SG::DBManager;
using SG::FEMSolver::AssemblyEigen;
// using SG::FEMSolver::CoordinateProcessor;
using SG::FEMSolver::ModelInfoQuery;
using SG::FEMSolver::Element::Hexa8SRIElementCalculator;
using SG::FEMSolver::Element::SolidMaterialParser;
using SG::FEMSolver::Element::Hexa8NodeShpFunc;

void Hexa8SRIElementCalculator::Initialize (SG::DataStructure::FEM::ElementTypes type)
{
    m_type             = type;
    auto integPointSet = m_integCalc.GetPoints ();
    m_shpFunc.Compute (integPointSet);
    m_shpFuncIncom.Compute (integPointSet);
}

void Hexa8SRIElementCalculator::CalcJ0Matrix (const SG::Algebra::Matrix& localcoord, _OUT SG::Algebra::Matrix& invJ0,_OUT Real& detJ0)
{
    // 计算J0, detJ0
    const auto tmp = 1.0/ 8.0;
    auto dNdxi_centor = Matrixd (3, 8, {
        -tmp, tmp, tmp, -tmp, -tmp, tmp, tmp, -tmp,
        -tmp, -tmp, tmp, tmp, -tmp, -tmp, tmp, tmp,
        -tmp, -tmp, -tmp, -tmp, tmp, tmp, tmp, tmp
    });
    const auto J0 = dNdxi_centor * localcoord.Transpose();
    invJ0 = SG::Algebra::inverse3X3 (J0, detJ0);
}

void Hexa8SRIElementCalculator::CalcBMatrix (const SG::Algebra::Matrix& localcoord, int iGauss, _OUT SG::Algebra::Matrix& B,_OUT Real& detJ)
{
    // 计算J
    const auto dNdxi        = m_shpFunc.GetShapeFunctionDerivative(iGauss);
    const auto J = dNdxi * localcoord.Transpose();
    // 计算invJ
    const auto invJ = SG::Algebra::inverse3X3 (J, detJ);
    // 计算dNdx
    const auto dNdx = invJ * dNdxi;
    // 利用dNdx计算B矩阵
    // 应变矩阵{epsx,epsy,epsz,gammaxy,gammaxz,gammayz}
    B = Matrixd (6, 24);
    for (std::size_t iNode{ 0 }; iNode < 8; ++iNode)
    {
        const auto i = iNode * 3;
        B (0, i)     = dNdx (0, iNode);
        B (1, i + 1) = dNdx (1, iNode);
        B (2, i + 2) = dNdx (2, iNode);
        B (3, i)     = dNdx (1, iNode);
        B (3, i + 1) = dNdx (0, iNode);
        B (4, i)     = dNdx (2, iNode);
        B (4, i + 2) = dNdx (0, iNode);
        B (5, i + 1) = dNdx (2, iNode);
        B (5, i + 2) = dNdx (1, iNode);
    }
}

void Hexa8SRIElementCalculator::CalcGMatrix (int iGauss,const SG::Algebra::Matrixd invJ0, Real J0OverJ, _OUT SG::Algebra::Matrix& G)
{
    // h1=1-xi^2 h2=1-eta^2 h3=1-zeta^2
    const auto dhdxi        = m_shpFuncIncom.GetShapeFunctionDerivate(iGauss);
    // 使用求逆犯规法计算计算dhdx
    const auto dhdx = invJ0 * dhdxi * J0OverJ;
    // 利用dhdx计算G矩阵
    G = Matrixd (6, 9);
    for (std::size_t iNode{ 0 }; iNode < 3; ++iNode)
    {
        const auto i = iNode * 3;
        G (0, i)     = dhdx (0, iNode);
        G (1, i + 1) = dhdx (1, iNode);
        G (2, i + 2) = dhdx (2, iNode);
        G (3, i)     = dhdx (1, iNode);
        G (3, i + 1) = dhdx (0, iNode);
        G (4, i)     = dhdx (2, iNode);
        G (4, i + 2) = dhdx (0, iNode);
        G (5, i + 1) = dhdx (2, iNode);
        G (5, i + 2) = dhdx (1, iNode);
    }
}

void Hexa8SRIElementCalculator::ComputeStiffness (const SG::Algebra::Matrixd& localCoord,
                                                  const SG::Algebra::Matrixd& D,
                                                  _OUT SG::Algebra::Matrixd& K)
{
    auto integPointSet = m_integCalc.GetPoints ();
    int nGauss             = integPointSet.size ();
    Matrixd B;
    Matrixd G;
    Matrixd Kuu (24, 24);
    Matrixd Kua (24, 9);
    Matrixd Kaa (9, 9);
    Matrixd invJ0;
    Real detJ;
    Real detJ0;
    CalcJ0Matrix (localCoord ,_OUT invJ0,_OUT detJ0);
    for (int iGauss = 0; iGauss < nGauss; ++iGauss)
    {
        CalcBMatrix (localCoord, iGauss,_OUT B,_OUT detJ);
        CalcGMatrix (iGauss, invJ0, detJ0 / detJ, _OUT G);
        Kua += B.Transpose() * D * G * detJ;
        Kaa += G.Transpose() * D * G * detJ;
        Kuu += B.Transpose() * D * B * detJ;
    }
    // 静凝聚消去内部自由度
    K = Kuu - Kua * Kaa.Inverse () * Kua.Transpose ();
}

void Hexa8SRIElementCalculator::ComputeInternalCoordTrans (const SG::Algebra::PointArray& pointSet, _OUT SG::Algebra::Matrixd& TS)
{

}
//单元组装
void Hexa8SRIElementCalculator::ComputeElementMatries (const SG::DBManager::DBServiceFactorySP dbServiceFactory)
{
    auto&       assembly       = AssemblyEigen::GetInstance ();
    static const DofSet dofSet{DOFDirections::DOF_Ux, DOFDirections::DOF_Uy, DOFDirections::DOF_Uz};

    CHEXA8Data element;
    std::vector<std::size_t> nodeIndexSet(8);
    std::vector<int> dofGlobalIndexSet(24);
    PointArray               pointSet (8);
    IdArray dispCoordIdSet(8);
    PSOLIDData property;

    Matrixd kg(24,24);

    //数据存储接口
    const auto& nodeService            = dbServiceFactory->get<INodeService> ();
    const auto& eleService             = dbServiceFactory->get<IElementService>();
    const auto& propertyService        = dbServiceFactory->get<IPSOLIDService> ();
    const auto& materialSummaryService = dbServiceFactory->get<IMaterialSummaryService> ();

    const auto  eleSize                = eleService->GetElementSize (m_type);
    for(std::size_t iEle{0}; iEle < eleSize; ++iEle)
    {
        //获取节点数据:节点索引、节点在Basic坐标系下的坐标、节点位移坐标系ID
        eleService->FindByIndex (iEle, _OUT element);
        ModelInfoQuery::QueryNode (nodeService,element.m_g, 8, pointSet, nodeIndexSet, dispCoordIdSet);
        //单元坐标系下节点坐标矩阵
        Matrixd coord (3, 8);
        for (std::size_t iNode{ 0 }; iNode < 8; ++iNode)
        {
            coord (0, iNode) = pointSet[iNode].m_x;
            coord (1, iNode) = pointSet[iNode].m_y;
            coord (2, iNode) = pointSet[iNode].m_z;
        }

        //获取单元自由度对应的全局自由度索引
        assembly.GetGlobalDOFId(nodeIndexSet, dofSet, _OUT dofGlobalIndexSet);

        //获取材料信息
        propertyService->FindById (element.m_pId, _OUT property);
        auto    matTypes = materialSummaryService->GetMaterialType (property.m_mId);
        //计算物理矩阵
        Matrixd TMat (3, 3);
        Matrixd D (6, 6);
        Real    rho;
        SolidMaterialParser::ComputeElasticMatrix (dbServiceFactory, property.m_mId, TMat, _OUT D, _OUT rho);
        ComputeStiffness (coord, D, _OUT kg);
        //组装刚度矩阵
        assembly.AccumulateStiffness (dofGlobalIndexSet, kg);
    }
}


// 计算位移输出坐标系下应力刚度矩阵
void Hexa8SRIElementCalculator::ComputeStressStiffness (const SG::DBManager::DBServiceFactorySP dbServiceFactory, int staticSubId)
{}


